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We consider joint probability distributions for the class of coupled Langevin equations introduced 
by Fogedby [H.C. Fogedby, Phys. Rev. E 50, 1657 (1994)]. We generalize well-known results for 
the single time probability distributions to the case of N-time joint probability distributions. It is 
shown that these probability distribution functions can be obtained by an integral transform from 
distributions of a Markovian process. The integral kernel obeys a partial differential equation with 
fractional time derivatives reflecting the non-Markovian character of the process. 
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I. INTRODUCTION 

In recent years, the connections between "continuous time random walk" (CTRW), which originated in the work of 
Montroll and Weiss [1 generalizing the idea of Brownian random walks, and fractional Fokker-Planck equations have 
been established. For a review we refer the reader to |2]. The solutions of these equations exhibit both super- and 
subdiffusive behaviour and are thus appropriate models for a large variety of transport processes in complex systems 
[3] . Recently, a connection between the velocity increment statistics of a Lagrangian tracer particle in fully developed 
turbulent flows and a type of CTRW has been introduced U . Here, a closure assumption on a hierarchy of joint 
velocity-position pdf 's derived from a statistical formulation of the Navier-Stokes equation leads to a generalization of 
Obukhov's random walk model [5] in terms of a continous time random walk. It allows for a successful parametrization 
of the single time probability distributions of velocity increments. However, there are different suggestions for the 
stochastic process of Lagrangian particles in turbulence, which are able to provide reasonable approximations for the 
single time velocity increment statistics. This example evidences that one has to introduce further quantities in order 
to distinguish between different stochastic models. 

For non-Markovian processes, the natural extension is the consideration of N-times joint probability distributions. 
It seems that for the class of CTRWs only single time probability distributions have been investigated so far. In that 
case fractional diffusion equations of the form 

~/(M) = Dl- a Lf(x,t) (1) 

can be derived. Here x denotes the random variable, L is a Fokker-Planck operator (for diffusion processes L — ) 
and oD\~ a is the Riemann-Liouville fractional differential operator (c.f. appendix A). The properties of this equation 
with regard to physical applications have been extensively discussed in the recent reviews [2], [BJ. In [7J Fogedby 
introduced a class of coupled Langevin equations, where he also considered a case which leads to an operator L 
including fractional derivatives with respect to the variable x, L — . A similar case has been studied by Meerschaert 
et al. |5], who made an extension to several dimensions introducing a multidimensional generalization of fractional 
diffusion, so-called operator Levy motion. This allows for a description of anomalous diffusion with direction dependent 
Hurst indices Hi defined by the relation < (xi(t) — Xi(t = 0)) 2 >~ t 2Hi . In jS] limit theorems of a class of continuous 
time random walks with infinite mean waiting times have been investigated. It is shown that the limit process obeys 
a fractional Cauchy problem. The emphasis again is put on single time distributions. 

The purpose of the present paper is to investigate multiple time probability distribution functions for the class 
of coupled Langevin equations introduced by Fogedby [7J, which have been considered to be a representation of a 
continuous time random walk. 

The paper is outlined as follows. In the next section we present the coupled Langevin equations considered 
by Fogedby [7] consisting of a usual Langevin process X(s) in a coordinate s and a Levy process representing a 
stochastic relation t(s). One is interested in the process X(t) = X(s^ 1 (t)). Fogedby .7] investigated the case where 
the processes X(s) and t(s) are statistically independent and showed how fractional diffusion equations of the form 
arise. Interesting results for the case where the processes are statistically dependent have been considered by 
Becker-Kern et al. [ID] leading to generalizations of the fractional diffusion equations ([!]). However, both publications 
are devoted to single time probability distributions. 

In section II we present a central formula, which relates the N-times probability distributions of X(t) to the pdf 's 
of X(s) via an integral transform, which is determined by the process t(s). In section III properties of the involved 
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Levy-stable process t(s) are considered leading to expressions for the pdf of the inverse process s(t). In section V we 
specify the moments for the case of a simple diffusion process. 



II. A CLASS OF NON-MARKOVIAN PROCESSES 



Starting point of our discussion is the set of coupled Langevin equations [7] for the motion of a Brownian particle 
in an external force field F in d=l dimensions (an extension to higher dimensions d > 1 is straightforward): 



dX(s) 
ds 

dt(s) 



F(X)+r,(s), (2) 

r(s) . (3) 
ds 

In this framework the random walk is parametrized in terms of the continuous path variable s, which may be 
considered eg. as arc length along the trajectory. X(s) and t(s) denote the position and time in physical space. The 
random variables r/(s) and r(s) are responsible for the stochastic character of the process. We are only considering the 
case of uncoupled jump lengths and waiting times such that r\ and r are statistically independent (coupled CTRWs 
have been considered in [TO]). The arc lenght is related to physical time t by the inverse function s = t~ l (t) = s(t). 
Thus, we have to assume r(s) > 0. We are interested in the process X(s(t)), i.e. the behaviour of the variable X as 
a function of physical time t. 

For the characterization of the process we introduce the two-times probability density functions (pdf) for the processes 

fi(x 2 ,s 2 ;x 1 ,s 1 ) = < 8{x 2 -X(s 2 ))S(x 1 - X( Sl )) > , (4) 
p(t2,s 2 ;t 1 ,s 1 ) = <S(t 2 -t(s 2 ))S(t 1 -t(s 1 ))> , (5) 
f(x 2 ,t 2 ;x 1 ,t 1 ) = <6(x 2 -X(s(t 2 )))5(x 1 -X(s(t 1 )))> . (6) 

Here the brackets < .. > denote a suitable average over stochastic realizations. For the sake of simplicity we restrict 
ourselves to n = 2. The generalization to multiple times is obvious. Both probability functions are determined by the 
statistics of the independent random variables r\ and t. 



A. The process X{s) 

We consider the case where ri(s) is the standard Langevin force, i.e. 77 is a Wiener process. In turn ^ becomes 
Markovian and fi(x 2 , s 2 ; xi, si) can be determined by solving the corresponding Fokker-Planck equation (FPE) for 
the conditional probability distribution P(x 2 , s 2 \ x\, s\): 

d ( d d 2 \ 

— P(x 2 ,s 2 I x 1 ,s 1 ) = i- — F(x) + —\P(x 2 ,s 2 \x 1 ,s 1 ) 

= L FP (x)P(x 2 ,s 2 I xi,si) . (7) 

The diffusion constant is set to 1 in the following. Due to the Markovian property of the process X(s) the joint pdf 
is obtained by multiplication with the single time pdf according to 

fi(x 2 , s 2 ; xi, 81) = P{x 2 ,s 2 I xi,si)/(xi,si) . (8) 

For a general treatment of the FPE we refer the reader to the monographs of Risken pTTJ and Gardiner [TJ] . 



B. The process t(s) 



The stochastic process t(s) is determined by the properties of r(s). The corresponding pdf's are denoted by p(t, s), 
p(^2, S2', ti, si). Furthermore, we shall consider r(s) to be a (one-sided) Levy-stable process of order a [7], |13j with 
< a < 1. As a result, the process t(s) is Markovian. Levy-stable processes of this kind induce the property of a 
diverging characteristic waiting time < t(s) > . Consequently the stochastic process in physical time t, given by the 
coupling of the Langevin equations |2]) and ([3| reveals subdiffusive behaviour. The specific form of p(t2, S2', t\, s\) will 
be given below. 

For a deeper discussion we refer to the review articles [B], [5] where the general relation between subdiffusive 
behaviour and diverging waiting times has been treated in detail. 
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FIG. 1: Sketch of the process t(s) which relates the arc length s to physical time t. Since the increment r(s) of eq.|3| is positive, 
the curve t(s) is monotonically increasing, implying the validity of the relation (14 1. . 



C. The process X(t) = X(s(i)) 

We are interested in the properties of the variable X with respect to the physical time t. Therefore, we have to 
consider the inverse of the stochastic process t = t(s): 



s = t~ 1 (t) = s(t) 

The stochastic process X(s(t)) then is described by the joint probability distribution 

f(x 2 ,t 2 ;x 1 ,t 1 ) =< 6(x 2 - X(s 2 ))S(s 2 - s(t 2 ))6(xi - X( Sl ))6( Sl - s(ti)) > 



(9) 



(10) 



The N-point distributions are determined in a similar way. Introducing the probability distribution h for the inverse 
process s(t), 



h(s,t) = <6(s- s(t)) > 
h(s 2l t 2 ;si,t 1 ) = < 8(s 2 - s{t 2 ))5{si - s{t{)) > 



(11) 



we can calculate the pdf of the process X(t) — X(s(t)) as a function of the physical time by eliminating the path 
variables s,-: 



/■OO /"OO 

f(x 2 ,t 2 ;x 1 ,t 1 ) = / cLsi / ds 2 h(s 2 ,t 2 ;s 1 ,t 1 )f 1 (x 2 ,s 2 ;x 1 ,s 1 ) 
Jo Jo 



(12) 



This relationship is due to the fact that the processes X(s) and t(s) are statistically independent. In that case, the 
expectation values in (10 1 factorize. Equation (12) can be generalized to N times. In fact, one may turn over to a 



path integral representation: 

f(x(t)) = J ■Ds(t)h(s(t))f 1 (x(s(t))) 
However, we do not investigate this path integral further. 



(13) 



The probability distribution h can be determined with the help of the cumulative distribution function of 
s(t). Since the process t(s) has the property (for s > 0) s 2 > si — > t(s 2 ) > t(si), one has the relationship 



e(s- s(t)) = i - e(t-t(s)) 



(14) 



Here, we have introduced the Heaviside step function: 0(x) = 1 for x > and Q(x) = for x < 0, Q(x = 0) = 1/2. 
The validity of eq.(14l becomes evident from an inspection of fig. 1: The function 0(s — s(t)) equals one in the 
region above the curve t — t(s), whereas Q(t — t(s) equals one in the region below the curve t = t(s). On the curve 

e(s-s(t)) = 1/2 = e(t-t( s )). 
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An immediate consequence is the following connection among the cumulative distribution functions of the processes 
t(s) and s(t): 

<e(« -*(<)) > = l- <Q(t-t( s )) > 
< e( S2 - s(t 2 ))e(«i - a(h) > = < (i - e(t 2 - t(s a )))(i - ©fa - t(si))) > 

= l- < e(fa - t(s 2 )) > - < e(h - t(si)) > 

+ < e(ta - t(s 2 ))Q(h - > • (15) 



Simple differentiation of eq.(15l yields the probability density function /i of the process s(t): 



h(s,t) = < Q(t-t{s)) > 

as 

h(s 2 ,t 2 ; Sl ,h) = ——<e(t 2 -t(s 2 ))e(t 1 -t(s 1 ))> . (16) 

Furthermore, since for t = 0we have the correspondence s = 0, the usual boundary conditions hold: 

h(s,0) = S(s) 
h(s 2 ,t 2 ;s 1 ,0) = h(s 2 ,t 2 )S(si) 
h(s 2 ,t 2 -> h;si,h) = S{s 2 - si)h(si,h) , (17) 



and can be verified from eq.(16l 



III. DETERMINATION OF THE PROBABILITY DISTRIBUTIONS p(s,t): LEVY-STABLE PROCESSES 

In the following we shall consider the joint multiple times pdf of the Levy-stable process ^ of order a. Simple 
integration of ([3| yields 

t( Si ) = I ' ds'T(s') , (18) 
Jo 

where we assume r(s) > 0. Additionally we consider the characteristic function for lu = iX. This defines the Laplace 
transform 



/>oo />oo 

Z(X 2 ,s 2 ;X 1 ,s 1 ) :=£{p{t 2 ,s 2 ;h,si)}= / dt 2 dhe- x ^- x ^p{ 

Jo Jo 



t 2 ,s 2 ;h,si) ■ (19) 



It will become clear below that working with Laplace transforms is more convenient for manipulating the pdf's of 
process ^ in the present context. 

A. One-sided Levy-stable processes: Single time 

At this point we have to introduce specific properties of the Levy-stable process. Levy distributions L a ^{x) are 
defined by two parameters |14J . |15j : a characterizes the asymptotic behaviour of the stable distribution for large x 
and hence the critical order of diverging moments. /3 characterizes the asymmetry. In the present case r > and 
the distribution is maximally asymmetric p(t < 0, s) = 0. This leads to (3 = 1. In the following we denote the Levy 
distribution L a> p(x) for (3 = 1 by L a (x). 

Let us motivate the consideration of Levy statistics. To this end we consider the characteristic function, which we 
write in the form: 

Z{\, a) = < er Xsl/ ° Jo ds '^'> > , (20) 
where a is a certain parameter. The choice Z(X,s) = Z(X a s) leads to a scale invariant pdf pit, s) = l/s 1 / a P(^7^) 
As a result, the characteristic function takes the form 



Z(X,s) = e- A ° s 



(21) 
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where we assume < a < 1. 

The probability distribution then becomes 

P(*>*) = ^M^) ' ( 22 ) 
where L a (i) denotes the one sided Levy stable distribution whose Laplace transform is £{L a (t)} — e _A °. 



B. Multiple times 

The joint pdf of the Levy process t(s) has been introduced in eq.([5|. Starting with this definition the derivation of 
the explicit expression for the pdf is straightforward and clearly reveals the Markovian character of this process. The 
characteristic function is given as Laplace transform of eq.([5]): 

/>oo />oo 

Z(X 2 ,s 2 ;X 1 ,s 1 ) = / dt 2 / dt ie ~ X2t2 ~ Xltl p(t 2 , s 2 ;t 1 ,s 1 ) 
Jo Jo 

= < e ~>*So 2d! "T{s')-\ 1 fi 1 ds'T(a') > ^ 

For further evaluating this expression we have to distinguish between the cases s 2 > s± and s± > s 2 . With a given 
ordering of s 2l s\ we can rearrange the integrals and write Z as a sum of two contributions: 

Z{X 2 ,s 2 ;Xi,si) = &{s 2 -s 1 )<e Js i Jo > 

. , -Ai / 32 ds'r(s') — (Ai+A 2 ) I " 2 ds'r(s') 

+6(si-s 2 )<e J =2 Jo > . (24) 

Here the expectation values factorize due to statistical independence of the increments r and can be expressed 
according to eq.(21 ): 

Z(X 2 ,s 2 ;X 1 ,s l ) = e( S2 - Sl )e' s ^ +x ^e~^- s ^ 

+6( Sl - S2 ) e -^(Ai + A 2 )" e -( Sl - S2 )A° _ ^ 

This is the characteristic function of the Levy process for multiple times. The appearance of the exponents (Ai + X 2 ) a 
is characteristic in this context and carries over to the pdf of the inverse process. We obtain the pdf p(s 2 ,t 2 ; Si,£i) 
after performing the inverse Laplace transform of eq. ( 25 1 . The result is 

■ (26) 

This expression explicitly exhibits the Markovian nature of the process. The conditional pdf p(t 2l s 2 \ti, Si) for s 2 > si 
is just: 

p(h, s 2 \h, Sl ) = 1 ~^L a ( h ~ t ' a ) ■ (27) 

(s 2 -Sl) 1/a \(S2-Sl) i/ "/ 

We remind the reader that L a (x) — for negative values of x. The expression for the joint pdf for multiple points is 
obvious. 



IV. THE PROBABILITY DISTRIBUTIONS h(s, t) 



The pdf's h(s,t), h(s 2 , t 2 ; si, t\) of the inverse process s — s(t) can be obtained from the pdf's of the process 
t = t(s) with the help of relationship eq.(16l. We shall consider the single- and multiple-time cases separately. Again, 
due to the simple form of the Levy distributions in Laplace space, we perform most of the calculations with Laplace 
transforms. 
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A. Single time 



Using the notation h(s,X) = £{h(s,t)} for the Laplace transform of h(s,t) with respect to t, the relation eq.(16l 
reads: 



Ms, A) 



. L < I e -AtW > = _ A I 

ds X ds X 



Z(s,X) 



The derivative with respect to s is easily performed with eq.(21 1 and leads to the solution h(s, A): 

h{s,X) = A a, - 1 e- sA ° 



(28) 



(29) 



This expression has already been derived in |7J — however without giving a 'simple physical argument'. Here the 
derivation is clearly based on eq.(14) which relates the Levy-stable process and its inverse. 
The inverse Laplace transform of eq.(29) is known and has been calculated in [T5] : 



Moreover, in |T7j the single time distribution h(s, t) has been identified as the Mittag-Lefiler distribution: 



ir , A (st a ) n 

h(s,t) = V A- '— 

i r(l + na 



(30) 



(31) 



Here we have obtained the pdf of s(t) for single times. Therefore, a complete characterization of the inverse process 
is given in this case. 

However in order to derive an evolution equation for the pdf of the process X(s(t)) we require an equation which 
determines h(s , t). 

From eq.(29) it is evident that h(s, A) obeys the differential equation 



J^( s > A) = X a h(s, A) 



(32) 



with the initial condition h(0, A) = A" 1 for s = 0. Hence, Laplace inversion yields a fractional evolution equation for 
h(s,t): 



-oD t 



l-a 



ds 



h(s,t) 



(33) 



The operator oD t ~ a denotes the Riemann-Liouville fractional differential operator, a possible generalization of integer 
order differentiation and integration to fractional orders (see Appendix B). For a discussion of fractional derivatives 
we refer the reader to [TBI. 



B. Multiple times 

The statistical characterization of the process s(t) for multiple times has been investigated from a mathematical 
point of view in the work of Bingham |17j already in 1971. He derived the following relationships for the moments 
< s(t N )...s(h) >: 

qn 1 

W . t, < s(t N )...8(h) > = 5yvv[*l(*2 " h)...(t N - tw-l)]"" 1 (34) 

oti...otN 1 (a) 



This equation can be obtained from the previous relation (16 1, which inferos the following relationship between the 
probability densities pit, s) and h(s, t): 

d . d . 

d 2 d 2 

-h(s2,t2;s 1 ,t 2 ) = t. — w—p(t2,s 2 ;t 1 ,s 1 ) 



q, i-iJ- V Z ' Z ' A ' ^ d Q 

OTlOT2 OS 2 OSl 

d N d N 

— -—h(s N ,t N ;...;si,t 2 ) = ^— p(t;v, sjv; si) • (35) 

Oti...Ot N OS N ...OSi 



7 



In the following we shall derive explicit expressions for these moments and show that instead of (34 1 fractional 
equations can be used for their determination. Based on eq.(16l and eq.(25l the derivation of an expression for the 
Laplace transform h(s 2 , A 2 ; s±, Ai) := £{h(s2, t 2 ; si, is obtained in a way analogous to the single-time case. 



We start by considering eq.(16l in Laplace-space: 



ft(*2,A 2 ;«i,Ai) = < — e -A**(«0jL e -Ai*(«i) > 



_d__d_ 1 
dsi ds 2 A 2 
d d 1 
dsi ds2 AiA 2 



1 

A? 



Using eq.(25) we can perform the derivatives of ■Z , (A 2 , S2', Ai, s\) with respect to s%, S2'- 

h(s 2 ,X 2 ;s 1 ,X 1 ) = S(s 2 - Sl ) A " " (A \ + , A2r + Ag e-^+^ a 

+6(S2 _ Si) (Ag)((Ai + A 2 )° - A?) e _ (Al+A2) . sle _ A g (s2 _ Sl) 
A1A2 

+Q(si _ S2 ) ^ A?) ^ Al + ^ ~ W e -(A 1+ A 2 )" S2e -A ? ( Sl - S;i ) 

Ai A2 



(36) 



(37) 



As a result we have obtained the Laplace transform of the joint pdf h(s2,t 2 ', s±,t±). Unfortunately, a closed form 
of the inverse Laplace transform could not be calculated. The given solution h can be readily used however to derive 
meaningful expressions which characterize the inverse process s(t). 



1. Moments of the inverse process 



In order to obtain further information about the process s(t) for multiple times we calculate the moments of the 
pdf. Let us first demonstrate how this can be achieved for the simple case < s(ti)s(t2) >■ This moment is defined 
from the pdf h(s2, i 2 ; s\, t\) as: 



< s(*i)s(t 2 ) > 



/>oo />oo 

/ dsi / ds 2 sis 2 h(s 2 ,t 2 ; si,*i) 
Jo Jo 

f poo pOO 

Cr 1 \ I dsi / ds 2 sis 2 h(s 2 , A 2 ;si,Ai) 



(38) 



where the last step follows by interchanging inverse Laplace transform and integration. The integrations with respect 
to si, S2 can be simply performed with the help of expression eq.(36). The result is: 



ds\ 



ds 2 sis 2 h(s 2 , A 2 ; si, Ai) = (Ai + A 2 ) 



\— a,— 1 \— a— 1 
-a J 1 , A 2 



A2 Ai 

Now the inverse Laplace transform leads to an analytical solution for < s(ti)s(t2) > (see Appendix B): 



(39) 



< s(tx)s(t 2 ) > = e(t 2 -ti) 
+e(ti - 1 2 ) 



—r — ri? Q + , 1 r- n t?tS F ( a, -a; a + 1: — 

r(2a+l) 1 T(a+ l) 2 1 2 V ' t 2 



1 



• 2a 1 4-Ot+a 771 1 1 -1 ^2 

. r 9 + =77 ttt tit, r a, —a a + 1 — 

r(2a+l) 2 r(a + l) 2 1 2 V *i 



(40) 



Here F(a, b; c; z) denotes the hypergeometric function (see e.g. Ch.15 in [19]). 



One notices that in the limit t 2 — ¥ ti expression (40) agrees with the second moment < s{t) 2 >: 

2 



< S (i) 2 > = C- 1 {J™ s a A 0, - 1 e- A '"dfl| 



2<v 



r(2a + l) 



where eq.(29l has been used. The simple single time moment < s(t) > is given as < s(t) >= C 1 {A a x } 



r(a+i) 



(41) 
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The calculation of higher order moments essentially follows the same steps. 

Furthermore, we introduce the operator (^^- + m ^ ne sense of the single-time Riemann-Liouville fractional 

differential operator: C{(^- + ^-^j g(ti,t 2 )} = (\\ + \ 2 )~ a g(\i, A2) (see Appendix A). An explicit expression in 
terms of an integral reads: 

Wi+W2 ) 9(t^) = m l ^"W*i -*',*,-*') . (42) 

Using this fractional differential operator, we are in the position to write down a simple recursion relation for 
arbitrary moments of h({si,ti\). The second moment eq.(39) reads: 

/ \ ~ a 

< s(h)s{t 2 ) > = ( — + — ) {< s(h) > + < s(* 2 ) >} . (43) 



0*1 0*2 , 

This immediately leads to (we assume t 2 > *i): 

< s(* 2 )s(*i) > = [ D- a {<s(t 2 - ti + £1) > + <*(*0>}]| 1=tl ■ (44) 
The explicit expression allows one to obtain the fusion rule 

lim < s(* 2 )s(*i) > = < s(h) 2 >= f 1 dt' t"*- 1 < s{h -*')> = 2 D7 a s{t 1 ). (45) 

r(a) 7 

The calculation of the third order moment < s(ti)s(t 2 )s(t 3 ) > along the same lines yields the result: 

< S(* 1 )S(* 2 ) S (*3) > ^ijt 1 + W 2 + W^) { < S ^)^ 2 ) > + < S (*l)<*3) > 

+ < s{t 2 )s{h) >} . (46) 

The third moment is obtained via fractional integration of the sum of second order moments. In the general case, the 
n-th order moment is calculated by fractional integration with respect to n times of the sum of all permutations of 
71—I order moments. 

Due to the representation of the fractional operator 

1 <?(*i, *2,* 3 ) = ^rr / dt' t /a - 1 g(t 1 -t',t 2 -t , ,t 3 -t / ), (47) 



.0*1 0*2 0W T(a) J 

we can derive the fusion rule 

1 Z"' 1 

lim < s(* 3 )s(* 2 )s(*i) > = — - / d*'*'"" 1 {< s(*i -*>(*i -*') > +2 < s(* 2 - *')s(*i - *') >} 
t 3 ^ti+o r(a) J 

= o^ Q {<s(*i)s(*i)>+2< S (*2-*i + *i)s(^i)>}t 1 =t 1 • (48) 

The fusion t 2 — > *i leads to 

< s(*!) 3 > = 3 D- a < s (<!) 2 > = 6D- a D- a < s(*0 > = 6 0J D t f " < sfa) > . (49) 
The n-th order generalization reads: 

< s(f) n > = n\ D; {n ~ 1)a < s(t) > . (50) 

This equation can also be derived directly from /i(s, A). Thus one can obtain a complete characterization of the 
process s(t) based on eq.(37l or eq.(36l respectively. Below, we shall show how to obtain these results on the basis of 
an evolution equation for the multipoint pdf h(si,t±; sn^n). 
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2. The structure of the N-time pdf 



From eq.( 16 1 one can derive the general form of the pdf h of the inverse process s(t) . The two times pdf reads (here 
we assume the case s 2 > si for simplicity) 

d d f* 1 f t2 
h(s 2l t 2 ;si,ti) = -7—-^— I dt[ / dt 2 p(t 2 - t x , s 2 - si) p(t x , si) 



dsi ds 2 Jq 

d r* 



dt' 1 h(s 2 -s 1 ,t 2 -t' 1 )p(t' 1 ,s 1 ) . (51) 



We define 



d f 1 

H(s 2 - si,t 2 - h;si ~ s ,h -t ) = -—— / dt' 1 h(s 2 - s 1 ,t 2 -t' 1 )p(t' 1 -t ,s 1 - s Q ). (52) 



d d 

h(s 3 ,t 3 ;s 2 ,t 2 ;s 1 ,t 1 ) = ^—^—j dt[l dt' 2 h(s 3 - s 2 ,t 3 - t' 2 ) p(t' 2 - s 2 - Si) p^, s x ) 



The form of the three times pdf is obtained in the same way and reads for s 3 > s 2 > si 

d d r 
dsi ds 2 Jo 



with a straightforward extension to the general case. 

With the help of eq.(52l this expression can be represented according to 



(53) 



d r t% 

K s a,ta;s 2 ,t 2 ;s 1 ,t 1 ) = — — / dt t H(s 3 - s 2 ,t 3 - t 2 ; s 2 - si,t 2 - p^, s x ) . (54) 

osi Jo 

Recursively, we may define higher order functions 

H N (sjv — Sff-\,tN — tjv-ij ti — to, s i — s o) 

= / dt' 1 H N ~ 1 (s N -s N -i,t N -t N -i;...;s 2 -si,t 2 -t' 1 )p(t' 1 -to,si,s ). (55) 

OS! Jo 

The integrals cannot simply be evaluated and the relations are formal. However, they show the underlying mathe- 
matical structure of the statistical description of the inverse process s(t). 



3. Fractional evolution equation 

In analogy to the single time case, where we have specified a fractional differential equation for h(s,t) 7 we now 
establish an evolution equation for h(s 2 ,t 2 ; si,tx). 
From eq.(37l it is evident that the following equation holds: 



d d_ 
dsi ds 2 



h(s 2 ,X 2 ;s 1 ,X 1 ) = -(Ai + X 2 ) a h(s 2 ,X 2 ;s 1 ,X 1 ) 



(56) 



with initial conditions 



ft(0,A 2 ;0,Ai) 



A" — (Ai + A 2 ) a + X 2 
AiA 2 



2 S2 



M»,A 2 ;0,A 1 ) . 

AiA 2 



(57) 



A common way to solve first order partial differential equations is the method of characteristics. Applying this method 
to eq.(56) with the given initial condition for each case , one obtains the correct expressions eq.(37l. Therefore eq.(56l 
determines the pdf in Laplace space. 
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Consequently, upon performing the inverse Laplace transform, we derive that h(s2, si, ti) obeys the fractional 
evolution equation 

J^r + J has been defined according to ( + j -F(i 2 , ii) : = 



(sfi - + zf^) (mT ^) F{t2,t\). The appearance of fractional time derivatives in eq.(58l reveals the non- 
Markovian character of the stochastic process s(t) and as a consequence of the coupled process - ? (s(t)). 



The extension of the above result to n times is straightforward: 



Again we want to emphasize that this single evolution equation with the proper initial condition sufficiently assets 
the pdf for multiple times. 

The above equation may also be used to calculate the moments < s(t/v)...s(ii) >, which already have been specified 



above. The fractional evolution equation (59) inferos the following relationship among the moments < s(ijv)...s(ti) >: 
( N d \ ( N d \ 1 ~° > 

\J2 of. J < s (*Jv)-s(ti) > = f ^2 gf J {< s(tjv-i)-s(*i) > +Permut} . (60) 



These equations are equivalent to the chain of equations (46) obtained by a direct inspection of the pdf's 



V. TWO-TIME MOMENTS OF THE DIFFUSION PROCESS 

In this last section we focus on the usual diffusion process, i.e. we consider the Fokker-Planck operator 

d 2 

L =w ■ < 61 > 

In this case, the moments are polynomials in s and we may directly use the results of the preceding session: 

< x(s 2 )x(si) > = 6(s 2 ~ + 6(si - s 2 )s 2 ■ (62) 
The corresponding moment with respect to time t is given by 

/>oo />oo 

< x(t 2 )x(ti) > = / / dsids 2 h(s2, t 2 \ si, ii) < x(s 2 )x(si) > . (63) 
Jo Jo 

The integrations can be performed by inserting the pdf h in Laplace space: 

C{< x(t 2 )x( tl ) >} = (Al + A2) " T ds s e-^rs = 1 . (64) 

AiA 2 J Q (Ai + A 2 j"AiA 2 

The inverse transform leads to the result 

<x(t 2 )x(t 1 )> = \ {Q(t 2 - t x )ff + 0fa - f a )ff } 
1 (a + 1) 

= e(t 2 - h) < s (h) > +e(h - 1 2 ) < s (t 2 ) > . (65) 

Similarly, we may calculate the moment < x{t 2 ) 2 x(ti) 2 >: 

< x(s 2 ) 2 x( Sl ) 2 > = s 2 si +26(s 2 - Sl )s 2 + 26( Sl - s 2 )s 2 2 . (66) 

This yields 

< x(t 2 ) 2 x(t!) 2 > = < s(i a )a(ti) > +29(t 2 - it) < s(ii) 2 > +20(ii - t 2 ) < s(t 2 ) 2 > . (67) 
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For the evaluation of < x(s 2 ) 2m x(si) 2n > we may use the properties of the moments of Gaussian processes which 
read for n > m: 

< x(s 2 ) 2m x( Sl ) 2n > = As™4 + Be(s 2 - Sl )si l - m 4 n + se(a x - s 2 )s^- m s^ . (68) 

The coefficients A, B, C can be evaluated by an application of Wick's theorem for Gaussian processes. 
The corresponding expression for the process X(t) becomes accordingly: 

< x(t 2 ) 2m x{h) 2n > = A < s^sih)" > +BO(t 2 - h) < s{h) n - m s(t 2 ) m > 

+B6(ii - 1 2 ) < s(t 2 ) n - m s(h) m > . (69) 

The calculation of the expectation values < s(t 2 ) 2m s(ti) 2n > has been discussed above. 



VI. CONCLUSION 



Up to now the discussion of continuous time random walks and the corresponding fractional kinetic equations has 
been focused on single time probability distributions only. On the basis of this pdf scaling behaviour of moments 
have been compared with experiments. However, more information has to be used in order to assign a definite 
stochastic process to a non-Markovian process. To this end we have considered multiple times pdf for a certain class 
of stochastic processes. 

Our approach is based on the framework of coupled Langevin equations ([2|,(|3| devised by Fogedby as a realization 
of a continuous time random walk. Here, the solution for the N-times pdf 's are given as an integral transform of the 
pdf 's of an accompanying Markovian process. We have shown that the non-Markovian character of this process can 
be traced back to the properties of the inverse Levy-stable process. 

The next step would be to compare these theoretical predictions with the behaviour of physical systems which reveal 
subdiffusive behaviour. To our knowledge multiple time statistics of such systems have not yet been investigated 
experimentally. This would be of considerable interest. We may expect that in some cases the consideration of 
multiple time statistics may lead to a more precise characterization of the underlying stochastic process. 

It is well-known, that for the single time case a fractional diffusion equation can be derived, which determines the 
pdf /(*, t), 

f(x,t)= / dah{8,t)fi(x,s) , (70) 
Jo 

as a solution of 

~f(x,t)= Dl- a L FP f(x,t) . (71) 

We would like to mention that a similar equation can be derived for the multiple times pdf f(x 2 , t 2 ]X\,ti). This will 
be discussed in a future publication. The present article is a starting point for the investigation of multiple times 
pdf 's of the coupled Langevin equations of Fogedby. 
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APPENDIX A: FRACTIONAL DIFFERENTIAL OPERATOR 



The Riemann-Liouville fractional integral is defined as a generalization of the Cauchy formula to real orders a: 

D7 a g(t) := -L [ - dt' 
4 y{ > T(a) J (t-t') 1 -" 

' t"- 1 *g(t) . (Al) 



T(a) 
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Here * denotes a Laplace convolution. Consequently performing the Laplace transformation is straightforward and 
yields the well-known result: 



£{oA~ Q 5(i)} = A- Q §(A) 



;From eq.( Al I the Riemann-Liouville fractional differential operator is obtained by simple partial derivation: 

D$- a g(t) := | Q D^g{t) . 



(A2) 



(A3) 



The extension of the fractional differential operator to two times ii,t 2 is now obtained in a way analogous to the 
steps above. 

First we define the fractional integral operator of two times in Laplace-space: 



d_ d_ 
dh + dT 2 



Furthermore the following equation holds: 



g{t u t 2 )\:= (Ai + A 2 )- a 5(Ai,A 2 ) 



(A4) 



POO POO i POO i 

/ dtj dt 2 e -Aiti-fcfa "tr'sih - h) = / dh e-^^ —tr 1 

Jo Jo T ( a ) Jo r (") 



= (A! + A 2 )- Q 



(A5) 



In physical time the fractional integral operator can thus be considered as an expression containing a two-fold 
Laplace convolution with respect to t\ and t 2 , denoted with **: 



d_ d_ 



l 



T(a) 



r(«) Jo 



dt[ / dt' 2 t?- l 8(t' 2 - t\) g(t 2 - t' 2 ,h - t[) 



(A6) 



Here we can distinguish between the cases t 2 < t\ and t 2 > t\ which results in eq.(47l The fractional differential 
operator of two times is then corresponding to eq.(A3l: 



d d 



g(h,t 2 ) 



dt^dt 2 ) \W 1 + df 2 



dU dt 2 , 

In the general N-times case the fractional integral operator takes the form of an N-fold convolution 

/ N 



d_ 



V 2 

\i=l 

with Laplace-transform 
C 



j(ti, ...,t N ) 



1 



T(a) 



t°T 1 S(t N - t N ^)...S(t 2 - fc) * ... * g(h, t N ) 



(A7) 



(A8) 



' N 

E 



dti 



g{tx, ...,t N ) 



/ n y a 

(E A M 5(Ai,-.-,Aat) 



(A9) 



APPENDIX B: CALCULATION OF MOMENTS 



Using the results of the previous section we can explicitly write the second order moment eq.(43) as convolution 
integrals: 



< s(ti)s(t 2 ) > = 



1 



r O) Jo Jo 



dt\ / dt' 2 tT~'5{t' 2 - t{) { 



1 



r(a + l) 



i 



r(a + l) 



(t 2 -tX} ■ (Bl) 



If we distinguish between the cases ti > t\ and t\ > t 2 in order to perform the integrations, we obtain: 

<«<*)-(*)> - ^ - { f(2o~+T)^ a + + 1) r dt> ^ - tT } 

+e(<1 - t2) { fX^TT) ^ + r (a) rU) C " ^ ~ °°} ' 

The integrals can be performed with Maple and lead to the hypergcomctric function F(a, b; c; z): 

/ 1 dt' t'^ih - t') a = -t*q F{a, ~a; a + 1; -) 
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